High-energy plasmon spectroscopy of freestanding multilayer graphene 
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We present several applications of the layered electron gas model to electron energy loss spectroscopy of 
free-standing films consisting of N graphene layers in a scanning transmission electron microscope. Using a 
two-fluid model for the single-layer polarizability, we discuss the evolution of high-energy plasmon spectra with 
N, and find good agreement with the recent experimental data for both multi-layer graphene with TV < 10, and 
thick slabs of graphite. Such applications of these analytical models help shed light on several features observed 
in the plasmon spectra of those structures, including the role of plasmon dispersion, dynamic interference in 
excitations of various plasmon eigenmodes, as well as the relevance of the bulk plasmon bands, rather than 
surface plasmons, in classifying the plasmon peaks. 
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I. INTRODUCTION 

Studying the relation between the plasmon spectra in thick 
graphite samples and those in carbon nanostructures with one 
or several layers has recently come to focus in carbon re- 
search. In that context, multilayer graphene (MLG) offers 
a particularly suitable system because of the precision that 
can be achieved in determining the number of carbon lay- 
ers N, as shown in the recent experimental study of plas- 
mon excitations in free-standing MLG by electron energy loss 
spectroscopy (EELS) in scanning transmission electron mi- 
croscope (STEM). 1,2 Those authors found that the so-called 
low-loss EEL spectra are dominated by two broad peaks that 
may be related to the ir and a + it plasmons of graphite, 
with the peak positions that change as the number N of lay- 
ers in the MLG increases. 1 ' 2 Similarly, plasmon spectra of 
graphene were also discussed in the context of recent exper- 
imental studies of thick layers of highly oriented pyrolytic 
graphite (HOPG) by ultrafast electron microscopy (UEM) in 
TEM 3 ' 4 and by inelastic x-ray scattering (IXS)r^ as well as 
in the experiments on multi-walled carbon nanotubes (MWC- 
NTs) using both EELS and IXS. 7 Given that in some of those 
studies the observed dependence of plasmon frequencies on 
the number of layers has been tentatively described in terms 
of the familiar concept of surface and bulk plasmons, it is de- 
sirable to adopt an analytically tractable theoretical model that 
can tackle this issue in an explicit and transparent manner. 

We demonstrate that a dielectric-response approach based 
on the layered electron gas (LEG) model, where each layer 
of carbon atoms is viewed as a two-dimensional electron 
gas (2DEG) of zero thickness, 8 - 9 provides useful insight in 
the effect of the number of layers on the plasmon spec- 
tra of MLG. The applicability of the LEG model to study- 
ing the high-frequency plasmon modes in MLG by means 
of EELS is justified for typical experimental conditions in 
STEM, where the target excitation is dominated by the in- 
plane transfer of momentum q of a fast incident electron 
that traverses a sequence of graphene layers, 1 ^* while elec- 



tron hopping between those layers may be neglected for exci- 
tation energies that exceed the interlayer coupling of t± « 
0.4 eV.— The LEG model has a long history in the study 
of plasmon excitations in various structures, including semi- 
conductor superlattices-iSii 2 . graphite intercalated compounds 
(GICs), 13 high-T c superconductors^ and MWCNTs^ 
A similar model has also been recently discussed in Ref. 6 
as a promising theoretical tool for extracting spectroscopic 
information on single-layer graphene (SLG) from their IXS 
data for thick graphite samples. Given the increasing number 
of recent developments in high-energy plasmon spectroscopy 
of graphene-based structures, 1-7 one is led to a conclusion 
that revisiting the old LEG model is well warranted. Due to 
its physical transparency and analytical flexibility, the LEG 
model offers new insight into various experimental observa- 
tions that complements the insight based on the results of 
computationally intensive ab initio approaches. 

The authors of Ref. 2 provided a theoretical discussion of 
their data by means of ab initio calculations of the loss func- 
tion for MLG in the optical limit, q = 0, based on the work 
of Marinopoulos et a/., 18 where HOPG was modeled by us- 
ing a supercell of carbon layers separated by an equilibrium 
distance of d w 3.35 A. In Ref. 2 an isolated SLG and iso- 
lated MLGs with N = 2 and 3 layers, having an interlayer 
separation d, were simulated by using supercells where those 
structures were periodically repeated with a separation be- 
tween them taken to be a multiple of d (typically fivefold). 
As a consequence of using supercells, the treatment of in- 
terlayer Coulomb interaction in ab initio calculations for the 
thus modeled SLG and MLGs is necessarily approximate, giv- 
ing rise to plasmon peaks in the loss function whose intensi- 
ties and the "precise peak positions depend on the separation" 
that was adopted within a supercell, as noted by the authors 
of Ref. 2 . Nevertheless, the ab initio calculations gave good 
qualitative agreement with the experimental EEL spectra for 
a SLG and MLGs with N = 2 and 3 layers. 2 Moreover, that 
agreement was taken as an indication that the optical limit of 
a loss function suffices for modeling the EEL spectra, even 
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though it was noted that those spectra were recorded so that 
the wave vector "q has a considerable in-plane component." 2 
Namely, while the 100 keV electrons traverse the MLG tar- 
gets undergoing negligible momentum transfer, the relatively 
large collection semiangle of 19 mrad in STEM ensured that 
the EEL spectra were recorded as being integrated over wave 
numbers up to q c ss 3.2 A -1 . 2 Furthermore, the authors of 
Refi 2 - did not consider any effects of the incident electron tra- 
jectory, even though such effects may give rise to a significant 
interference in plasmon excitations at different carbon layers, 
which is critically dependent on proper use of the interlayer 
separation. Finally, we mention that elements of phenomeno- 
logical modeling of the EEL spectra are introduced in the ab 
initio calculations by their use of a spectral broadening of 1.5 
eV to smooth out the resulting theoretical spectra, and hence 
improve comparison with the experiment* 2 - 

A possibly interesting alternative theoretical discussion of 
the experimental EEL spectra 2 could be based on a contin- 
uum dielectric model, which was used with notable success 
for multilayer carbon nanostructures, different from the MLG, 
by treating them as slabs of finite thickness described by a 
frequency-dependent dielectric tensor in the optical limit^^r— 
However, as pointed in Ref. 2 , most applications of such an 
anisotropic dielectric slab (ADS) model were restricted to the 
EELS of curved nanostructures with nonpenetrating, or aloof 
electron trajectories in STEM, so that excitations of the bulk 
plasmon modes were suppressed with respect to the surface 
modes. 19 22 It should be mentioned that the ADS model was 
used by Crawford 23 to study the stopping and deflection of 
swift ions passing through the bulk of HOPG, but it is true that 
full application of the ADS model to plasmon spectroscopy of 
MLG with penetrating electron trajectories is yet to be under- 
taken. On the other hand, while the ADS model is perceived 
as suitable for carbon nanostructures with a large number of 
layers, difficulties arise if one attempts to reach the limit of a 
one-atom-thick layer by starting from a dielectric slab of fi- 
nite thickness. 22 Namely, it was found that the slab thickness 
strongly affects the coupling of surface plasmons at the oppo- 
site sides of the slab, 24 and hence the authors remarked that 
achieving the limit of a single-layer is not straightforward in 
the ADS model because an "arbitrary choice of the effective 
dielectric thickness" has to be made.— 

It is expected that the LEG model can provide useful addi- 
tional insight into the issues discussed in the above two para- 
graphs regarding (a) the interlayer Coulomb interactions, (b) 
the effects of large in-plane momentum transfer, (c) the dy- 
namic interference effect due to electron trajectory, (d) the 
lack of available applications of the ADS modeU£~— for pene- 
trating electron trajectories, and (e) the single-layer limit. An- 
other benefit coming from the LEG model is revealed in the 
limit of an infinite periodic lattice of identical layers, giving 
rise to a particularly transparent description of the bulk plas- 
mon bands that may be of interest for studying the plasmon 
spectra of HOPGi 10 ' 11 ' 13 Finally, it is particularly convenient 
that the LEG model involves separate treatments of the inter- 
layer Coulomb interactions and the intralayer dynamics by 
using a wave vector-dependent, non-interacting polarizabil- 
ity function Xo(<l>k') of SLG as an independent input quan- 



tity. Obviously, whether fine details in the EEL spectra can be 
successfully modeled depends on the availability of a good- 
quality single-layer polarizability xo(<T w )> but the effects of 
the increasing number of layers within MLG are expected to 
be robustly reproduced by the LEG model. 

Over the past few years, sophisticated models have been 
developed for xoCt w ) that are appropriate for describing the 
low-energy excitations (up to, say, 1-2 eV) 25 in the vicin- 
ity of the K points in the Brillouin zone of SLG, where the 
conduction and valence tt electron bands may be approxi- 
mated by Dirac cones with zero gap.^— The intra-band tt 
plasmon excitations (sometimes called sheet plasmons) were 
probed in this energy range by means of high-resolution re- 
flection EELS (HREELS) of epitaxial graphene under heavy 
doping conditions* 2 ^ 2 ^ 2 - giving dispersion relations that 
were studied theoretically by means of Xo(q, u>), which in- 
cluded the effects of plasmon damping 33 and plasmon-phonon 
coupling! 2 ^ 4 - However, plasmon excitations at such low ener- 
gies are hardly accessible in the TEM experiments because 
of the presence of the zero loss peak that masks the spectral 
features up to about 2 eV. 35 On the other hand, high-energy 
plasmon spectra, which were observed in thick graphite sam- 
ples by means of EELS, may be associated with the inter- 
band excitations of both the it and a electron bands with the 
gaps of about 4 eV and 14 eV, respectively (see Ref3 and 
the references therein). As far as the modeling of such high- 
energy spectra in SLG is concerned, we note that the few re- 
cent improvements of ab initio calculations of the optical di- 
electric function of graphene did not quite agree with regard 
to the role of the excitonic effects associated with the a elec- 
tron bands. 36-38 Therefore, there is some value in the trans- 
parency offered by a phenomenological two-fluid model for 
high-energy excitations of the tt and a electrons in carbon. 39 
Accordingly, we adopt here such a model for SLG, while be- 
ing aware that subtle features due to low-energy, intra-band tt 
electron excitations are inaccessible in the EELS anyway.— A 
tradeoff to using a phenomenological xo(<T w ) as an input to 
the LEG model is that the resulting analytical tractability may 
help reveal how the main features in the EEL spectra of MLG, 
i.e., the tt and a + tt plasmon peaks, evolve as the number of 
layers increases from N = 1 in SLG to N — >• oo in HOPG. 
In particular, the LEG model will enable us to analyze the 
role played in the spectra by the formation of the bulk plas- 
mon bands in HOPG, and to show that the concept of surface 
plasmons has limited applicability in the present context. 

After presenting theoretical details of the LEG model in 
the following section, we discuss the comparison of our cal- 
culations with several experiments using a phenomenological 
model for the single-layer polarizability, which is followed by 
our concluding remarks. Note that we use Gaussian electro- 
static units, and we set H = 1. 



II. THE LAYERED ELECTRON GAS MODEL 

In a typical STEM experiment operating at the voltage of 
100 kV, 2 the momentum transfer of the incident electron is 
close to zero, so we shall use a straight-line trajectory while 
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neglecting relativistic effects, 35 We use a Cartesian coordinate 
system {r, z} with r = {x, y} and assume that graphene lay- 
ers occupy planes z n = (n—Tjd, where n = 1,2, ... ,N and d 
is the interlayer spacing. The induced potential in the system, 
<&i n d(r, z, t), may be expressed via its Fourier transform with 
respect to the in-plane coordinates (r — > q) and time (t — > cj) 
as 



where Pn(cj) is the probability density for loosing the energy 
cj, given by 



(Ze) 2 j _ 
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$ ind (q,z,cj) = — cr„(q,w)e ?|z z "l, 
^ — i y 



(i) 



where a n (q, w) is the Fourier transform of the induced charge 
density (per unit area) on the nth layer, which may be written 
in a self-consistent field approximation as 



<?„(q,w) = -e 2 xo(q,w) $cxt(q, z n ,u) + $ in d(q, z n , w) 



,(2) 



with $ ext (r, 2,t) being the external potential. From the 
charge density of the incident electron, p cxt (r,z,t) — 
ZeS(r— vut) 5(z— v±t), where Z = — 1 and V|| and v± 
are the velocity components parallel and perpendicular to the 
graphene planes, respectively, we find 



$cxt(q, z,u) 



4wZev\ 



(qv x ) + (w-q-V||)" 



a i(w-q-V||)a/i 



(3) 



Thus, by using Eqs. ([TJ and (0) in Eq. (12), we obtain a matrix 
equation for a n (q, cj), 



= -ZeV{q)x(q,u)K(q,uj-q ■ V\\)ip n (oj-q • vy), 
where 

x„„, ( q , w) = <W + (i - <W) fte)x(q,w)e~ 8,l|n '~ n| ,(5) 



^i{n — l)ujdfv± 



(6) 



(7) 



and x(q,w) = Xo(q,w)/e(q,a;), where e(q,w) = 1 + 
V (q)xo{<\, cj) is the intra-layer dielectric function of SLG 
with V(q) = 2ire 2 jq being the Fourier transformed Coulomb 
interaction in two dimensions (2D). 

In the next step we solve Eq. (@]i for <7 n (q, cj) by inverting 
the matrix A4 with elements given in Eq. (0. Substituting this 
solution into Eq. (Q3 enables one to express $i n d(r, z, t) via an 
inverse Fourier transform, so that the total energy lost by the 
incident electron may be evaluated from^ii 

oo oo 

E loss = - J dt J d 2 r J dzp cxt (r, z,t) — $ ind (r, z, t).(8) 



Using symmetry properties of the function xo(q 5 w ) at zero 
temperature, one may write Eq. ([H) as 



oo 

-Eloss = J ducjP N {cj), 



(9) 



^J^K 2 (q,cj-q- n )^[V(q) X (ci,cj)Q N (q,cj)}, 



(10) 



Q N (q,uj) = 

En=i E£=i <(«-q - v,|) (-M" 1 )^, *M"-q ■ V||) 



(11) 



Note that the density Pn(uj) will be directly compared with 
the experimental EEL spectra of MLG with finite N, which 
were taken under the normal electron incidence^ So, setting 
V|| = in Eqs. ( [TUt and ( fTTT i and invoking the near-isotropy 
of graphene's polarizabilityr^ xo(q, = Xo(q, w), renders 
the angular integration in Eq. ( TTOb trivial. The remaining in- 
tegration over the wave numbers should go up to q c w 3.2 
A -1 ^ but we found that no difference occurs in the final re- 
sults for Pn(u>) if the upper limit is extended to oo because 
the kinematic factor K 2 (q,uj) in Eq. (TTOb is strongly peaked 
at q = w/v± <C q c for frequencies of interest here, cf. Eq. ([6j. 

Further note that the factor [Ai ) , in Eq. (fTTT > is a (q, cj) 
dependent element of a matrix Mr 1 that is inverse to the ma- 
trix M. defined in Eq. (0. Thus, the quadratic form Qn{q., cj) 
in Eq. ( fTTT l can be relatively easily obtained using symbolic 
computation software for, say, N < 10. For example, whereas 
for a single-layer Qi(q, cj) — 1, we obtain from Eq. ( fTTT l for a 
bilayer graphene 

n ( x l-F(g)x(g^)e~ gd cosM/^) n ~ 

g2(9 ' w) - 2 — r^w&K^ — ' ( } 

While the analytical results for Qn become increasingly 
cumbersome with increasing N, we note that a relatively sim- 
ple expression for Pjv may be obtained in the limit N — » oo, 
which we shall denote by P m , corresponding to an electron 
traversing a sufficiently thick slab of HOPG, such that the end 
effects may be neglected if Ndq 1. In that case, the matrix 
equation in Eq. (0 may be solved by using Fourier series with 
a wavenumber k in the direction perpendicular to the graphene 
planes, giving 



PooH 



1+5 g 



V r (g)xa(g,w) 



(13) 



where 



5(g,A:) 



sinh^d) 



cosh(gc?) — cos(fcd) 



(14) 



is the Coulomb structure factor for an infinite periodic lattice 
of layers. 9 Note that the density P^cj) will be directly com- 
pared with the experimental EEL spectra of HOPG, which 
were taken under the normal electron incidence, 3 so we again 
set vii = in Eq. ( fT3b . 
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Finally, in order to identify the exact nature of plasmons 
giving the most prominent contributions to the spectral den- 
sities Pat(w) or P 00 (w), it is worthwhile analyzing the eigen- 
modes of the underlying MLG, which are obtained by set- 
ting the damping rates in xofe w) to zero. Thus, for plasmon , — , 
modes in the case of finite N, one has to solve the equation *u 
e N (q, oj) det {M.) = 0, giving frequencies at which the factor 
3 [V(q)x(q, w)Qjv(q, w)] in Eq. (IToT > becomes singular. On 3" 
the other hand, eigenmodes in an infinite periodic lattice of Oh 
2DEG layers are obtained by solving the equation 



l + S(q,k)V(q) X o(q,U J )=0, 



(15) 



with k as a parameter. 10, 1 1 By letting < k < tt jd and using a 
one-fluid model for xo (<Z, <*>), one obtains a band of dispersion 
relations for the so-called bulk plasmon modes that propagate 
with the wave numbers q parallel to the layers. 9-12 Similarly, 
by using a two-fluid model for Xo(<li w ) in Eq- (EK one ob- 
tains two such bands for < k < tt jd that correspond to 
the tt and a + tt plasmons in the bulk of HOPG. In the case 
of a semi-infinite lattice of equally spaced identical layers of 
2DEG, it was shown that a plasmon mode may arise with the 
dispersion relation outside the plasmon band(s), only if there 
is a mismatch between the background dielectric constants in 
the lattice and the nearby space»i2iii Hence, this kind of sur- 
face plasmon, which is localized near the boundary layer of a 
semi-infinite lattice, is not expected to exist in the N — > 00 
limit of the LEG model for HOPG placed in vacuum or air. 



III. RESULTS AND DISCUSSION 

A three-dimensional (3D) version of the phenomenological 
two-fluid polarizability function was used in the ADS model 
to build a dielectric tensor in the optical limit with suitable 
Drude-Lorentz parameters^ for modeling of the EEL spec- 
tra of multilayer fullerene molecules 20 and MWCNTsj^LZ 2 . In 
other applications, such a 3D version of the two-fluid model 
was used to study the variable degree of the sp 2 hybridiza- 
tion for applications in different carbon materials 41 and the 
in-plane plasmons in HOPG,— as well as to deduce the op- 
tical conductivity of graphene in order to calculate Casimir 
forces between graphene layers. 43 However, we need here a 
strictly 2D version of the two-fluid model with suitable Drude- 
Lorentz parameters, similar to that used to describe plasmon 
excitations in single-layer fullerene molecules 44,45 and single- 
wall carbon nanotubes (SWCNTs)4& One way of deriving 
such a polarization function for SLG could proceed from a 
2D, two-fluid hydrodynamic model with Thomas-Fermi and 
Dirac (TFD) interactions, 40 ' 47 which enabled a semiquanti- 
tative comparison with the plasmon dispersion relations that 
were observed in the EEL spectra of SWCNTs over a broad 
range of wavelengths in the axial direction.— 

We note that clear distinction should be made between the 
3D and 2D Drude-Lorentz models in the sense that the former 
class of models usually employs only frequency-dependent di- 
electric functions, whereas the latter class necessarily involves 
nonlocal effects, i.e., a dependence on the in-plane wave num- 
ber that reflects incomplete Coulomb screening by an electron 
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FIG. 1: (Color online.) Probability density Pn(lj) (in 1/eV) ver- 
sus energy loss ui (in eV), evaluated from Eq. d 1 0b for N = 1, 2, 
5, and 13 graphene layers [smooth solid (yellow, green, blue, and 
pink) curves], along with the corresponding experimental EEL spec- 
tra from Ref.~ [noisy (gray) curves]. 



gas in 2D.— A formal connection between the two types of 
models may be established by considering a thin film, where 
both symmetric and antisymmetric coupling of surface plas- 
mons at the opposite sides of the film arise, in addition to the 
bulk plasmon modes. Keeping in mind that, in realistic appli- 
cations to one-atom thick layers, only the surface density, n, 
of electrons may be defined unambiguously, taking the zero- 
thickness limit of a film leaves the lower-energy, symmetric, 
in-plane surface plasmon as the only observable excitation 
mode, characterized by a typical square-root plasmon disper- 
sion of the form to ~ y/nq in a quasifree 2DEG. 24 Whether 
such a plasmon mode should be referred to as a surface plas- 
mon, or simply an intrinsic plasmon mode of a 2DEG is a 
matter of semantics when it comes to one-atom-thick layers. 
Hence we adopt from Ref4°- a planar version of the 2D 

two-fluid hydrodynamic model that gives xo — Xt* 
for SLG, where 



(0) 



X i °Hl^) = 



nlq 2 /ml 



s 2 q 2 + ui 2 r -lo(lo + i-y v ) ' 



(16) 



with n°, m*, s v , u) vr , and 7„ being the equilibrium surface 
number density of electrons, effective electron mass, acoustic 
speed, restoring frequency, and the damping rate in the vi\\ 
fluid (where v = a, tt), respectively. Note that the restoring 
frequencies for the in-plane electron excitations are related to 
the tt — > tt* and a — > a* inter-band transitions, 44 which were 
found to dominate the in-plane loss function of SLG in the 
optical limit at energies close to 4 and 14 eV, respectively. 18 
On the other hand, the terms involving the acoustic speeds in 
Eq. ( fToT i arise from the TFD interactions in the hydrodynamic 
model, 40 but their contribution to the EEL spectra turns out 
to be negligible in the present context because the kinematic 
factor K 2 (q,uj) in Eqs. ( fTQb and ( fT3T > is strongly peaked at 
q = oj/v± and because s v <C v±. 

Taking to* to be the free electron mass and using the un- 
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perturbed surface electron densities of SLG, n° « 38 nm -2 
and n® = 3n^. sw 115 nm -2 , we treat the remaining parame- 
ters in Eq. (TToT l as adjustable. The best fit to the experimental 
EEL spectra 2 is found for bj^ r = 4.08 eV, io ar = 13.06 eV, 
7^ = 2.45 eV, and 7 CT = 2.72 eV. The results for P n {uj) 
were thus computed from Eq. ( TTOb with 2V = 1, 2, 5 and 
13, and are compared in Fig[T]with the experimental curves 
from Fig. 1(e) of Ref. 2 , corresponding to 1, 2, 5, and > 10 
graphene layers (We found that our curve TV = 13 provides 
the best fit with their curve labeled ">10L".). We note that the 
experimental curves were taken under the same acquisition 
conditions, thereby enabling a direct quantitative comparison 
among them. 2 On the other hand, apart from adjusting the ar- 
bitrary unit of the experimental curves to the absolute unit of 
our Pn(u>), no relative scaling took place among the theoreti- 
cal curves with different numbers of layers. One notes that the 
experimental curves are well reproduced by the LEG model 
for energies u> > 3 eV and N < 10, both in magnitude and 
in the shape of spectra. Regarding the discrepancy at uj < 3 
eV, we note that the present version of the two-fluid model in 
Eq. ([Tol l, along with the neglect of interlayer tunneling, is only 
expected to work in MLG for high-energy excitations, but also 
that the complete vanishing of the experimental spectra at en- 
ergies below 3 - 4 eV may be a consequence of the method 
used in their subtraction of the zero-loss peak. 2 

Furthermore, while the above values of the parameters u„ r , 
w CTr , 7tt, and 7 CT are quite close to those listed in the Table 
2 of Ref. 19 for the in-plane dielectric function of a 3D, two- 
fluid model of ADS, features seen in the spectra in FigQ]are 
relatively robustly reproduced by theoretical curves for other 
choices of these parameters. In particular, the damping rates 
are strongly affected by the presence of impurities or defects, 
which serve as scattering centers for charge carriers in individ- 
ual carbon layers. While for graphene on a substrate the con- 
centration of impurities strongly varies from sample to sam- 
ple, the freestanding graphene may be relatively clean in the 
case of one or few layers, but thicker samples may contain 
increased amounts of impurities and defects. Thus, assign- 
ing smaller values, or possibly allowing for frequency- de- 
pendent damping rates 7^ and 7^ as in Ref. 49 , could improve 
the agreement for energies around 10 eV in Fig.Q] where the 
experiment exhibits an almost complete depletion of the spec- 
tra. On the other hand, the width of the high-energy peak at 
about 27 eV for N > 10 is not well reproduced by the present 
choice of parameters, but agreement may be improved if one 
allows for higher damping rates due to increased density of 
impurities. 2 This point will be taken up later in the discussion 
of Fig. 

The most important trends seen in Fig.Q]are that the tt plas- 
mon peak position moves from about 5 eV to about 7 eV as 
N increases without significant changes in its shape, whereas 
the a + tt plasmon peak at about 15 eV for N = 1 evolves 
through the development of a plateau between 15 and 27 eV 
for N = 5, to be dominated by a peak at about 27 eV for 
N = 13, with a growth in magnitude that exceeds the growth 
of the tt plasmon peak. The lowest peak positions that oc- 
cur for N = 1 and the highest peak positions that occur for 
N > 10 in Fig.Q~]were associated in Ref. 2 with the surface and 
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FIG. 2: (Color online.) Probability density Pn(lj) (in 1/eV) versus 
energy loss to (in eV), evaluated from Eq. J 1 0b for N = 1 [lower (red) 
curves] and 2 [upper (blue) curves] graphene layers with non-zero 
damping (smooth thin solid lines) and vanishing damping (dashed 
lines), along with the corresponding experimental data [noisy (grey) 
curves] and the ab initio curves (smooth thick solid lines) from Ref. 2 . 



bulk plasmon frequencies of graphite, respectively, for both 
the 7r and a + tt plasmons. Similarly, we shall see in Fig. [6] 
that the a + tt plasmon contributions at about 15 and 27 eV in 
the EEL spectra of HOPG were also associated in Ref.— with 
the surface and bulk plasmons, respectively. In the following, 
we shall use the LEG model to discuss how adequate those 
associations are. 

In Fig. [2] we analyze the details of the experimental EEL 
spectra of Eberlein et al. for N = 1 and 2, corresponding to 
a SLG and a bilayer graphene (BLG), which are reproduced 
from Fig. Q] (noisy curves) along with the theoretical curves 
from the LEG model for P\(uj) and P2(co) from Eq. ( TlOb 
(thin solid lines) and with two theoretical curves, taken from 
Fig. 3(b) of Ref. 2 (thick sold lines), which are based on their 
ab initio calculations. The most prominent feature of the ex- 
perimental curves is the asymmetry in the a + tt plasmon 
peaks, which is well reproduced by the LEG model, but not 
by the ab initio calculations by Eberlein et al. 2 (We note, how- 
ever, that such asymmetry was observed in the ab initio calcu- 
lations based on a GW and Bethe-Salpeter equation approach 
by Trevisanutto et al£ty We emphasize that the long tails in 
those peaks, which extend at energies beyond 15 eV, are a con- 
sequence of the integration over q in Eq. ( fTOb that captures the 
effect of dispersion of the a + tt plasmon, which is neglected 
when optical limit of the loss function is used. To further elu- 
cidate this point, we have recalculated the curves for Pi (ui) 
and P2(cj) from Eq. <TT0b using the same parameters of the 
two-fluid model in Eq. ( fT6l ), except that the damping rates of 
both a and tt electrons were made very small. The resulting 
curves are displayed in Fig. [2] by the dashed lines, showing 
a steep raise at the restoring frequencies u) ur , along with the 
tails that extend for u> > w vr , giving a very good approxima- 
tion to the experimental curves at those frequencies. 
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This observation may be easily rationalized for N = 1 by 
referring to the two-fluid model for Xoili u ) m Eq. (O in the 
limit 7 CT = 7^ — > 0, which was discussed in the Appendix C 
to Ref. 40 . When used in Eq. (TlOb for Pi(w), this limit con- 
verts the factor 9 [V(q)x(q, w)] into a superposition of two 5 
functions, S(u 2 — iu±(q)), where to±(q) are the eigenfrequen- 
cies describing hybridization between the a and tt electron 
fluids in a SLG, which are obtained by solving the equation 
e(g, oj) = as 

u2± = ^+ny + r$ ± ^ ^+^i-ni y + nlQli (17) 
where we define 

w„(g) = V^r + s^ 2 (18) 
and 

n„(g) = v /2^e 2 nO (Z /m* (19) 

for the i^th fluid. 40 The eigenfrequencies w+ (g) and ui- (q) are 
labeled as the a + tt and the tt plasmon modes of SLG, and the 
corresponding dispersion relations are shown in Fig.[3]by the 
upper and lower dotted lines, respectively. For the purpose of 
the present discussion, it suffices to expand the expressions in 
Eq. ( PTTb to the first order in q, so that i>j\ k, uj 2 r + tt 2 (q) and 
w 2 w uj 2 r + f2^(g), describing the long wavelength behavior 
of plasmon frequencies of the non-interacting a and tt fluids in 
SLG. Using these expressions for w 2 . (g) in the delta functions 
S (uj 2 — oj± (q) ) makes the integration over q in Eq. ( TT~Qb trivial, 
giving two contributions to Pi(w) of the form ~ n1/(u 2 — 
w 2 ,,.), which essentially capture the behavior of the tails seen 
inFig.|2]forw 2 - w 2 r > 2ire 2 n%ujl (m> ± ). 

Furthermore, it is interesting to comment on the simi- 
larity between the experimental spectra for BLG and SLG, 
seen in Fig. [2] In particular, their high-frequency tails at 
uj > 20 eV seem to imply a scaling ratio of about 2:1, 
as remarked by Eberlein et air This is best rationalized 
by resorting to the eigenmode analysis for BLG by setting 
7 ff = — > in Eq. ( fTOb for P2(w). Then, it can be 
shown that the factor 3 [V(q)x(q, w)Q 2 (g, w)] becomes sin- 
gular at four frequencies, which are given by the expressions 
in Eq. (ITTb when frequencies il„(q) in Eq. (fT9l are replaced 
by (q) = fl v (q) a/1 ± exp(— qd). The corresponding dis- 
persion curves are shown in Fig.|3]by the dashed lines, two in 
the upper {a + tt plasmon) group and two in the lower {tt 
plasmon) group. However, the weakly (quadratically) dis- 
persing lower curves in each plasmon group contribute to 
P 2 (uj) with negligible spectral weights because of the small- 
ness of the factor ujd/v± at the frequencies of interest here. 
Namely, by setting cos (uid/v±) ss 1 in Eq. (1121 1. one obtains 
Q 2 (q,uj) ~ 2/ [1 + V(q)x(q, w) exp(— qd)], so that the dom- 
inant spectral weight in P2(w) comes mostly from the upper 
(linearly dispersing) curves in the ct+tt and tt plasmon groups, 
seen in Fig. [3] for BLG. This may be further simplified by 
noting that the kinematic factor K 2 (q,oj) is strongly peaked 
at q = u)/v±, so that exp(— qd) sa 1, and hence the fac- 
tor 3 [V(g)x(g, w)Q2(g, £*;)] in Eq. (|ToT > for P 2 (o;) becomes 
— 9f{l/ [1 + 2V(q)xo(q, w)]}, corresponding to a SLG with 
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FIG. 3: (Color online.) Dispersion curves for the n (lower group) and 
<t+7t (upper group) plasmons in MLG, obtained from the LEG model 
with N = 1 [dotted (red) lines] and N = 2 [dashed (blue) lines], 
as well as N = 5 [only the upper group is shown by (black) solid 
lines]. The upper and the lower edges are shown for both the n and 
<7 + 7r plasmon bands in HOPG [thick (orange) solid lines]. The (dark 
cyan) diamonds show the experimental plasmon peak positions in the 
EEL spectra of HOPG from Ref. 50 , whereas the (black) squares show 
the peak positions obtained by ab initio calculations for graphite in 



doubled a and it electron densities. Accordingly, the slopes 
of the upper dispersion curves for BLG in the a + tt and tt 
plasmon groups are seen in Fig.[3]to be about twice the slopes 
of the corresponding dispersion curves for SLG in the limit 
of long wavelengths. As a consequence, the peak positions in 
p2(w) occur at about the same frequencies as those in Pi(w), 
whereas the high-frequency tails in Pa(w) are approximated 
by ~ 2n° / (cj 2 — u) 2 r ), confirming their ratio of 2: 1 relative to 
the tails in Pi(u>). 

The above analysis of SLG and BLG shows that various 
eigenmodes of the underlying structure enter the EEL spectra 
with weights that are strongly dependent on the incident elec- 
tron speed. This dynamic aspect of the LEG model, which 
is not covered in the ab initio calculations^ is related to the 
interference effect in plasmon excitations at different carbon 
layers within a MLG, and it may be quantified by the factor 
ojTm that appears in the expression for Qjy(q, uj), cf. Eqs. ©, 
© and (flTT i. where Xjy = Nd/v± is the time it takes the 
incident electron to traverse the full thickness of MLG (ne- 
glecting relativistic effects). Therefore, the interference effect 
will be negligible for fast incident electrons in STEM at fre- 
quencies of interest in the low-loss EELS of MLGs with a 
few layers giving cjTjv -C 1. The smallness of the factor 
cjTn = Nduj /v± is equivalent to the limit d — » 0, giving a 
picture where all layers collapse onto a single sheet or, equiv- 
alently, the in-plane plasmon modes in all layers oscillate in 
phase, which explains the similarity among the spectra and 
the approximate scaling of their high-energy tails as Pjv oc N 
for N ~ 1, In other words, the EEL spectra of MLGs with 
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FIG. 4: (Color online.) Probability density Pm(oj) (in 1/eV) versus 
energy loss uj (in eV), evaluated from Eq. d 1 01 > for TV in the range 
from 1 to 15 with unit steps, using the same parameters as in Fig.Q] 
[solid (black) curves]. Also shown are the results of ab initio calcu- 
lations from Ref. 2 for graphene with N = 1,2, and 3 layers [dashed 
(red, blue, and green, respectively) curves]. 



N = 2,3,... are dominated by the peaks that are derived from 
the 7r and a + tt plasmons of SLG, and hence they should not 
be associated with the surface plasmons of graphite. 2 

When N > 10, the factor wT/v = Nduj/v± may no longer 
be small, and hence the interference effect may become sig- 
nificant in sufficiently thick MLGs. Obviously, the onset of 
the interference effect will occur at a smaller N if the inci- 
dent electron moves at a lower speed. Moreover, because of 
the difference in the relevant frequency ranges, one expects 
that the interference effect will be more prominent for the 
a + tt than for the tt plasmons. This may explain why the 
a + tt plasmon peak in Fig. Q] undergoes a more substantial 
change with increasing N than the tt plasmon peak. To il- 
lustrate this point, we show in Fig. [4] a detailed evolution of 
the EEL spectra -P/v(u;) from Eq. ( flOl l for N going from 1 
to 15 in unit steps, which were calculated with the same pa- 
rameters of the two-fluid model in Eq. ( [ToT l as those used in 
Fig- 03 (For the sa ke of further comparison with ab initio cal- 
culations, we also show in our Fig.|4]the results from Fig. 3(b) 
of Ref. 2 for the in-plane spectra of graphene with N = 1,2, 
and 3 layers.) One confirms in Fig.|4]that in the LEG model 
the 7r plasmon peak moves continuously from about 5 to about 
7 eV as N increases, whereas the spectral structure associated 
with the a + tt plasmon appears to evolve as a superposition 
of two peaks, one at about 15 eV and the other at about 27 
eV, which do not move much as N increases but rather have 
their weights strongly dependent on N. One may further see 
from Fig. [4] that the contribution from the peak at about 15 
eV dominates the spectra for N < 5, and is still visible for 
N < 10, whereas the spectra with N > 10 are dominated 
by the peak at about 27 eV. In particular, it appears in Fig. [4] 
that the weights of those two peaks are approximately equal 
for MLG with N — 6 layers, giving a value of the relevant 



interference factor of ujTq ss 0.5. Such a peculiar composi- 
tion of the a + tt plasmon peak for N ~ 5 - 6 deserves further 
analysis. 

One expects that the plasmon eigenfrequencies in MLG 
with finite N will fall into two disjoint groups, each giving N 
dispersion curves that correspond to the a + tt and tt plasmons 
of SLG. All the curves within each group become degener- 
ate as q — !• and they approach the corresponding restoring 
frequencies of SLG, uj ar and ui nr , but they spread out at fi- 
nite q, forming two quasibandsjiSiii In Fig. [3] we show such 
a quasi-band for N = 5 in the a + tt group (the correspond- 
ing quasi-band in the tt group is not shown to avoid cluttering 
of curves in that group). Similar quasibands may be seen in 
the case of MWCNTs with, e.g., N = 30 walls in Fig. 4(a) 
of Ref. 16 based on a single-fluid model for a electrons, and 
N = 10 walls in Fig. 1 of Ref. 17 based on a two-fluid model 
for both a and tt electrons. One sees in Fig. [3] that the four 
lower curves in the a + tt quasi-band for N — 5 are quadrati- 
cally dispersing, whereas the uppermost curve is linearly dis- 
persing with a slope that is about five times the slope of the 
corresponding N = 1 curve for small q. It can be shown that, 
by setting cos (uid/v±) ~ 1 as in the case with A" = 2, the 
uppermost curve in the a + tt quasi-band for N = 5 gives 
the dominant contribution to the spectrum Ps(u;). The large 
initial slope of that dispersion curve, along with its tendency 
to saturate for higher q values, gives rise to a high-energy tail 
in Ps(o;) that forms a bump at about 27 eV, which is clearly 
seen in the experimental spectrum for N = 5 in Fig. Q] As 
A" further increases, this saturation of the dispersion curves 
within a quasi-band becomes more prominent because of a 
well-defined upper bound, 16 which is best illustrated by re- 
sorting to the eigenmode analysis of an infinite periodic lattice 
ofSLGs. 

In the limit JV ^ oo we recover the case of HOPG, which 
exhibits two bands of dispersion relations corresponding to 
the a + tt and tt bulk plasmon modes. They are obtained by 
letting < k < tt jd and solving the equation in Eq. ( fTBT ) 
with Xo(q,w) from the two-fluid model of Eq. ([Tol l where 
la = 7?r —> 0*£~— The upper and the lower edges of 
those bands are defined by k — and k = tr/d, respec- 
tively, and the corresponding dispersion relations are given 
by uj — w± p (g) and u = wlg w (?) with "+" for the a + tt 
plasmon band and "-" for the tt plasmon band, where expres- 
sions in Eq. ( fTTl ) are to be used with the frequencies f2„(g) 
in Eq. (O replaced by fiy(q) = Q, v (q)y/coth(qd/2) and 
n l ° w (q) = n„(<?)v / tanh(gd/2). Those banded ges are shown 
by the four thick solid (orange) curves in Fig. [3] which all ex- 
hibit parabolic dispersions in the long wavelength limit. 

When ij -> 0, one can show that the lower band edges in 
Fig. [3] approach the restoring frequencies of SLG, wJP w (0) = 
uj ar and wl? w (0) = u> nr , whereas the upper band edges ap- 
proach the values W£ P (0) that are given by uj± in Eq. ( fTTl i 
when the frequencies f2„ (q) in Eq. ( fl9] i are replaced by Q vp = 

A:TTe 2 N2 /m*, corresponding to the bulk plasmon frequency 
in the z^th electron fluid with an effective volume density of 
N% = n° u /d. Using our parameters from Fig. Q] we find 
il^p ~ 21.7 eV and Vt^ p w 12.5 eV, which are quite close 
to the plasmon frequencies listed in the Table 2 of Ref^ for 
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the in-plane dielectric function of the ADS model. Hence, 
we obtain for the long wavelength limits of the upper edges 
of the bulk a + tt and tt plasmon bands in Fig. [3] the values 
w+ p (0) » 27.7 eV and w" p (0) 7 eV, respectively, which 
are remarkably close to the peak positions seen in the Pn(cj) 
curves in Figs. Q] and |4] for N > 10. This may be rational- 
ized by referring to the N ~ > oo limit of the LEG model, 
given in Eq. ([T3T l. where we see that the smallness of the fac- 
tor kd = ujd/v± in S(q, k), cf. Eq. ( [Pil l implies that the dom- 
inant spectral weights in P oc (o;) come from the upper edges 
(defined by k = 0) of the bulk a + tt and tt plasmon bands of 
HOPG. 

Hence, for finite q, one expects that the peak positions of 
both the 7r and a + tt plasmons in the EEL spectra of HOPG 
should closely follow the dispersion relations of the upper 
edges of the corresponding plasmon bands. This is corrob- 
orated in Fig. |3]by a comparison with the peak positions taken 
from Figs. 3 and 5 of Ref. 50 , which were deduced from the ex- 
perimental EEL spectra of HOPG. (We note that the overdis- 
persion of the a + tt upper band edge relative to the experi- 
mental points at larger q values in Fig. [3] may be reduced if a 
finite damping rate is introduced in the two-fluid model.) In 
addition, we display in Fig. [3] the peak positions obtained in 
RefJ£ from ab initio calculations for graphite, showing good 
agreement with both the experimental data 50 and the upper 
edges of the plasmon bands. Therefore, the previously men- 
tioned association of the plasmon peaks in the EEL spectra of 
thick samples of MLG or HOPG with the bulk plasmon modes 
of graphite 2,3 should be made more specific by stressing that 
it is the upper edges of the bulk plasmon bands that give the 
prevalent contributions to the EEL spectra of such structures. 

On the other hand, for incident electrons at lower speeds, 
one expects that increased contributions in the plasmon spec- 
tra of HOPG may also come from the portions of the bulk plas- 
mon bands with lower energies. In that context, we mention 
that Fig. 1 in Ref. 51 shows the results of a HREELS experi- 
ment on HOPG using a 200 eV incident electron beam, where 
contributions from the entire a + tt plasmon band are seen 
to strongly depend on the electron reflection angle. In partic- 
ular, a prominent peak was observed in the spectra at about 
13 - 14 eV that did not move much when the reflection angle 
was changed. 51 While that peak was associated by the author 
with a surface plasmon, its presence may be alternatively ex- 
plained as resulting from a van-Hove type of singularity in the 
plasmon density of states at the lower band edge of the a + tt 
plasmon band. 15 Clearly, further investigation into the features 
of high-energy spectra of HOPG in HREELS is warranted. 

It is worthwhile mentioning that the tt plasmon band was 
probed in a recent experiment by Hambach et a/., 5 who used 
IXS to measure the frequency dependent dynamic structure 
factor (DSF) of HOPG for momentum transfers beyond its 
first Brillouin zone. The experimental data from Fig. 3(a) 
in Ref. 5 are reproduced in our Fig. [5] along with the results 
from their ab initio calculations, showing the frequency de- 
pendence of a region near the tt plasmon peak for a range 
of perpendicular wave numbers, tt < kd < 3tt, with a fixed 
in-plane wave vector component of about 0.37 A -1 . In partic- 
ular, the peak position is seen to oscillate between the values 
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Energy loss [eV] 

FIG. 5: (Color online.) Dynamic structure factor [solid (blue) lines] 
is evaluated from Eq. J20b and appropriately scaled to compare it 
with the corresponding experimental IXS data (points) from Ref. 5 
with a fixed in-plane wavenumber of q = 0.37 A -1 for the indicated 
values of the perpendicular wavenumber k. Also shown are the re- 
sults of ab initio calculations from Ref. 5 [dashed (red) lines]. 



of about 6 and 7 eV, which are commensurate with the extent 
of the tt plasmon band in our Fig. [3] when one sets q = 0.37 
A -1 and invokes graphene's nearly isotropic in-plane polariz- 
ability. Furthermore, one can show that the DSF of an infinite 
periodic lattice is proportional to the factor 



V(q)xo(q,u) 



l + S(q,k)V(q) X o(q,u) 



(20) 



which appears in Eq. (TT3T > with S(q, k) defined in Eq. (fl4l) . 
By using the same parameters as those used in Fig. Q] one 
can easily evaluate from Eq. ( f20b a set of spectral curves for 
HOPG with q = 0.37 A" 1 for a range of kd values, which 
are seen in Fig. [5] to reproduce all the main features of the 
corresponding experimental data and ab initio results for DSF 
from Ref. 5 , including the tt plasmon peak position variation 
with kd, and its near disappearance for kd = 2tt. We take that 
this comparison demonstrates both the validity of the concept 
of plasmon bands in HOPG and the ability of the LEG model 
to describe them in a simple analytical manner. 

Regarding the applicability of the N — > oo limit of the LEG 
model in the context of EELS in STEM, we obviously require 
ujTn 3> 1 because, in view of the strong peaking of the kine- 
matic factor at q = lj/v± in P 00 (a;) from Eq. 1131 . this con- 
dition amounts to Ndq ^> 1, which guarantees that the end 
effects are negligible. When ujTn ~ 1, the interference is 
strong and the end effects in such MLG may not be neglected, 
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Energy [eV] 

FIG. 6: (Color online.) The EEL spectra versus energy loss (in eV), 
evaluated from Poo(w) in Eq. d 1 3fc [thick (black) solid line], along 
with the corresponding experimental spectrum of HOPG from Ref. 3 
[the (violet) chain curve], which is theoretically decomposed into a 
laser-generated peak [thin (black) solid line], n plasmon peak [thin 
(black) dash-dot line], a "surface" a + n plasmon peak [thin (black) 
dashed line], and a "bulk" a + n plasmon peak [thin (black) dotted 
line]. 



making it necessary to use Pat(w) from Eq. (fTOb . even when 
N ^> 1. To illustrate this point, we mention that, when at- 
tempting to model in Fig. Q] the experimental EEL spectra for 
the thickest MLG from Ref. 2 with "> 10" layers, we could 
not obtain nearly as good a fit with P^w) as we did with 
Pjv(w) for N = 13. The fact that the interference factor in 
this case takes the value U1T13 m 1 points to a need to look 
elsewhere for the experimental EEL spectra of thick enough 
MLG, where we may test the applicability of Poo(oj) from Eq. 
AD to HOPG. 

In that context, we find it instructive to model the recent 
experiment by Carbone et al.,— who used their UEM in the 
so-called static mode to obtain the EEL spectra of a slab of 
HOPG, having such a thickness that ujT n sa 20. In Fig. 
we compare the experimental curve from Fig. 2 of Refi 3 - for 
HOPG with our result for the (appropriately scaled) proba- 
bility density P oc (a;) from Eq. dT3l . where we used the two- 
fluid model for Xo(l> <*>) given in Eq. ( T3"6i >. In order to achieve 
the best fit with the experiment, we used the same parameters 
for the P QO (w) curve in Fig. [6] as those used for the Pat(cj) 
curves in Figs. [T] and @] except for 7^, which was increased 
to 7 CT = 9.52 eV to take into account that the damping rate 
for a electrons may be increased due to more abundant impu- 
rities or defects in a thick slab of HOPG than in a few-layer 
MLG. Our curve for Poo(oj) is seen in Fig. [6] to reproduce 
the experiment very well for energies > 10 eV due to the in- 
creased parameter 7^, whereas the tt plasmon peak is only 
qualitatively reproduced in shape, but its position fits the ex- 
periment well. The discrepancy between our model and the 
experiment at energies < 10 eV may be partially ascribed to 



the presence of both the zero-loss peak and the photogener- 
ated charge carrier plasma excitations at around 2.4 eV, which 
were not subtracted from the experimental curve in Ref.—. We 
also show in Fig.|4]a phenomenological decomposition of the 
experimental curve into four theoretical curves, which was 
performed in Ref A Those curves describe contributions from 
several processes that were associated by the authors with: the 
laser-generated region peaked at about 2.4 eV, the tt plasmon 
peaked at about 7 eV, the surface a + tt plasmon peaked at 
about 15 eV with a low weight, and the bulk a + tt plasmon 
peaked at about 27 eV with a much larger weight than its sur- 
face counterpart. 3 While such a decomposition of the a + tt 
plasmon contributions is corroborated by the above analysis, 
we believe that it would be more adequate to label the peaks 
at 15 and 27 eV as contributions related to the lower edge 
(which carries a signature of the SLG's a + tt plasmon mode) 
and the upper edge of the bulk a + tt plasmon band in HOPG, 
respectively. 



IV. CONCLUDING REMARKS 

We have shown that a simple and straightforward imple- 
mentation of the LEG model 91 1 to MLG gives analytical ex- 
pressions for the spectra of the energy lost by a fast electron 
under normal incidence that provide an explicit account of 
the effect of increasing number N of graphene layers on the 
high-energy plasmon excitations in MLG. Those expressions 
are given in terms of a single-layer polarizability Xo(<7>w) 
that may be modeled independently from the LEG model. 
By adopting a phenomenological, two-dimensional, two-fluid 
model for xo{q,w)> which includes Lorentz parameters suit- 
able for describing the dominant tt — > tt* and a —> a* inter- 
band transitions in SLG, we found good agreement with the 
plasmon spectra from four independent experiments. In par- 
ticular, the shape of such spectra was well reproduced for 
MLG with N < 10, where the most dramatic change oc- 
curs in the tt and a + tt plasmon peaks. 2 In addition, it was 
shown that the experimental EEL spectra for both MLG with 
N > 10 2 and a thick slab of HOPG 3 may be well reproduced 
by the same model by allowing for increased damping rates in 
Xo(?,w). 

Specifically, the LEG model reproduced the experimentally 
determined peak positions, 2 which were found to move from 
about 5 to about 7 eV for the tt plasmon and from about 15 to 
about 27 eV for the a + tt plasmon as the number of layers 
in MLG grows from N ~ 1 to N > 10. 2 By referring to the 
eigenmode analysis of the underlying MLG structures, both 
the plasmon peak positions and the experimentally observed 
similarity among the EEL spectra for N = 1,2,3,... were 
found to be derived from the plasmon spectra of SLG, based 
on the smallness of the factor ujd/v± under typical experi- 
mental conditions (low loss EELS and the 100 keV incident 
electrons with d rj 3.35 A -1 ). 2 Hence no association with 
the surface plasmons of graphite seems to be justified for the 
peaks at about 5 and about 15 eV for the tt and <t+tt plasmons, 
respectively, in MLG with few layers. In addition, the high- 
frequency tails in the experimental EEL spectra were exactly 
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reproduced by the LEG model for N = 1 and 2, and were 
identified as arising from the plasmon dispersion and from in- 
tegration over a large range of the in-plane wave numbers, 
commensurate with the experimental conditions. 2 

On the other hand, for sufficiently thick MLG, such that 
Ndtu/v± ^> 1, the limit of an infinite periodic lattice of SLGs 
showed that the peaks at about 7 and about 27 eV come pre- 
dominantly from the upper edges of the tt and a + it plasmon 
bands in the bulk of HOPG, respectively, again because of the 
smallness of the factor ud/v±. This seems to be corroborated 
by a comparison with the experimental dispersion relations 
found from the EEL spectra of graphite. 50 Moreover, it was 
shown that the N — > oo limit of the LEG model also gives 
good account of a dynamic structure factor that was measured 
by IXS of HOPG, 5 reproducing the movement of the tt plas- 
mon peak between the corresponding tt plasmon band edges. 

When Nduj/v± ~ 1, the a + tt plasmon peak is seen to 
evolve through a sequence of broad structures between 15 and 
27 eVr 2 which result from the development of a quasiband of 
plasmon dispersion curves for MLG with N ~ 10 that is a 
precursor to a fully developed bulk a + tt plasmon band of 
HOPG. Furthermore, having in mind that the underlying semi- 
infinite lattice of SLGs for N — > oo does not support surface 
states in HOPG placed in the air or vacuumrSiii one may con- 
clude that any remaining trace of the peak at about 15 eV for 
the a + tt plasmon in the EEL spectra of HOPG is likely to be 



related to the lower edge of its bulk a + tt plasmon band* 3 - and 
not to a surface plasmon of graphite. 

Given that small variations in the free parameters, used in 
the model adopted for xo(q> w )> do not corrupt the agreement 
found with the four independent experiments , 2 ' 3 ' 5 ' 50 we sug- 
gest that using the LEG model for high-energy plasmon ex- 
citations in MLG, HOPG, and other multilayer carbon nanos- 
tructures indeed presents a versatile and robust alternative to 
other theoretical approaches. The novel applications of these 
analytical models, presented in this work shed light on several 
aspects of the problem at hand, including (a) the role of plas- 
mon dispersion in the spectra integrated over the wave num- 
ber, (b) the role of the dynamical interference factor Ndu/v± 
in determining which eigenmodes of the underlying MLG 
structure have prevalent spectral weights, and (c) the relevance 
of the bulk plasmon bands, rather than surface plasmons, in 
classifying the observed plasmon peak frequencies. 
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